Food competition analysis

Author

Max Lindmark

Published

January 25, 2024

Load packages

home <- here::here()

library(INLA)
library(tidyverse)
library(tidylog)
library(RCurl)
library(sdmTMB)
library(RColorBrewer)
library(devtools)
library(patchwork)
library(ggstats)
library(ggh4x)
library(sdmTMBextra)

# Source map-plot
source_url("https://raw.githubusercontent.com/maxlindmark/cod-interactions/main/R/functions/map-plot.R")
#source(paste0(home, "/R/functions/map-plot.R"))

Read data & prepare data

d <- read_csv(paste0(home, "/data/clean/aggregated_stomach_data.csv"))

# Calculate relative prey weights (saduria and benthos)
d <- d %>% 
  drop_na(group) %>% 
  drop_na(oxy) %>% 
  rename(oxygen = oxy) %>% 
  mutate(tot_weight = rowSums(select(., ends_with('_tot'))),  
         benthic_weight = amphipoda_tot + bivalvia_tot + gadus_morhua_tot +
           gobiidae_tot + mysidae_tot + non_bio_tot + 
           other_crustacea_tot + other_tot + other_pisces_tot + platichthys_flesus_tot +
           polychaeta_tot + saduria_entomon_tot) %>% 
  rename(saduria_weight = saduria_entomon_tot,
         flounder_density = fle_kg_km2,
         large_cod_density = mcod_kg_km2,
         small_cod_density = scod_kg_km2) %>% 
  mutate(tot_rel_weight = tot_weight / (pred_weight_g - tot_weight), 
         benthic_rel_weight = benthic_weight / (pred_weight_g - tot_weight),
         saduria_rel_weight = saduria_weight / (pred_weight_g - tot_weight)) %>% 
  dplyr::select(-ends_with("_tot")) %>% 
  dplyr::select(-predator_latin_name, date) %>% 
  # scale variables
  mutate(fyear = as.factor(year),
         fquarter = as.factor(quarter),
         fhaul_id = as.factor(haul_id),
         depth_sc = as.numeric(scale(depth)),
         oxygen_sc = as.numeric(scale(oxygen)),
         density_saduria_sc = as.numeric(scale(density_saduria)),
         flounder_density_sc = as.numeric(scale(flounder_density)),
         large_cod_density_sc = as.numeric(scale(large_cod_density)),
         small_cod_density_sc = as.numeric(scale(small_cod_density)))

Quick explore

Sample size

d %>%
  group_by(species) %>% 
  summarise(n = n())
group_by: one grouping variable (species)
summarise: now 2 rows and 2 columns, ungrouped
# A tibble: 2 × 2
  species      n
  <chr>    <int>
1 Cod       5259
2 Flounder  3851
d %>%
  group_by(species, quarter) %>% 
  summarise(n = n())
group_by: 2 grouping variables (species, quarter)
summarise: now 4 rows and 3 columns, one group variable remaining (species)
# A tibble: 4 × 3
# Groups:   species [2]
  species  quarter     n
  <chr>      <dbl> <int>
1 Cod            1  3012
2 Cod            4  2247
3 Flounder       1  2081
4 Flounder       4  1770

Fit models

Groups are: small cod, large cod and flounder. Response variables are: saduria_rel_weight, benthic_rel_weight or total weight. The latter is only for adult cod, because essentially all prey are benthic for small cod and flounder.

# This is the reason we don't do total weight for flounder and small cod
d %>%
  filter(tot_rel_weight > 0) %>% 
  group_by(group) %>% 
  mutate(ben_prop = benthic_rel_weight / tot_rel_weight) %>% 
  summarise(mean_ben_prop = mean(ben_prop))
# A tibble: 3 × 2
  group     mean_ben_prop
  <chr>             <dbl>
1 flounder          0.978
2 large cod         0.592
3 small cod         0.956

Covariates are: ~ 0 + fyear + fquarter + depth_sc + spatial + spatiotemporal random fields + density covariates. For saduria, we use saduria also in interaction with cod and flounder. For cod we use small cod because large and small cod are very correlated. For benthic and total prey, we instead use oxygen, more as a proxy, as the interaction variable

Compare models with different spatial terms and evaluate AIC

  • Update, I now do 5-fold cross validation to select models, see script food_competition_crossvalidation.qmd

Main models

pred_flounder_sad <- list()
pred_flounder_ben <- list()
pred_cod_sad <- list()
pred_cod_ben <- list()
coef_sad <- list()
coef_ben <- list()
res_sad <- list()
res_ben <- list()
random_sad <- list()
random_ben <- list()
range_sad <- list()
range_ben <- list()

for(i in unique(d$group)) {
  
  dd <- filter(d, group == i)
  
  mesh <- make_mesh(dd,
                    xy_cols = c("X", "Y"),
                    cutoff = 10)

  ggplot() +
    inlabru::gg(mesh$mesh) +
    coord_fixed() +
    geom_point(aes(X, Y), data = dd, alpha = 0.2, size = 0.5) +
    annotate("text", -Inf, Inf, label = paste("n knots =", mesh$mesh$n), hjust = -0.3, vjust = 3) + 
    labs(x = "Easting (km)", y = "Northing (km)")
  
  ggsave(paste0(home, "/figures/supp/mesh_", i, ".pdf"), width = 17, height = 17, units = "cm")
  
  
  # Saduria model
  
  if(unique(dd$group) == "flounder") {
    
    m_sad <- sdmTMB(saduria_rel_weight ~ 0 + fyear + fquarter + depth_sc + oxygen_sc + 
                      small_cod_density_sc*density_saduria_sc + 
                      flounder_density_sc*density_saduria_sc,
                    data = dd,
                    mesh = mesh,
                    family = tweedie(),
                    spatiotemporal = "IID",
                    spatial = "off",
                    time = "year")
    print(i)
    sanity(m_sad)
    print(m_sad)
    
  } else {
    
    m_sad <- sdmTMB(saduria_rel_weight ~ 0 + fyear + fquarter + depth_sc + oxygen_sc + 
                      small_cod_density_sc*density_saduria_sc + 
                      flounder_density_sc*density_saduria_sc,
                    data = dd,
                    mesh = mesh,
                    family = tweedie(),
                    spatiotemporal = "IID", 
                    spatial = "on",
                    time = "year")
    print(i)
    sanity(m_sad)
    print(m_sad)
    
    }

  
  # Benthic model
  
    if(unique(dd$group) %in% c("large cod", "small cod")) {
      
      
      m_ben <- sdmTMB(benthic_rel_weight ~ 0 + fyear + fquarter + depth_sc + 
                        small_cod_density_sc*oxygen_sc + flounder_density_sc*oxygen_sc,
                      data = dd,
                      mesh = mesh,
                      family = tweedie(),
                      spatiotemporal = "IID", 
                      spatial = "off", # spatial off
                      time = "year")
      print(i)
      sanity(m_ben)
      print(m_ben)
      
    } else {
      
      m_ben <- sdmTMB(benthic_rel_weight ~ 0 + fyear + fquarter + depth_sc + 
                        small_cod_density_sc*oxygen_sc + flounder_density_sc*oxygen_sc,
                      data = dd,
                      mesh = mesh,
                      family = tweedie(),
                      spatiotemporal = "IID", 
                      spatial = "on",
                      time = "year")
      print(i)
      sanity(m_ben)
      print(m_ben)
      
    }
       
   
  # Spatial and spatiotemporal random effects
  d_haul <- dd %>%
    distinct(haul_id, .keep_all = TRUE)

  preds_sad <- predict(m_sad, newdata = d_haul)
  preds_ben <- predict(m_ben, newdata = d_haul)

  random_sad[[i]] <- preds_sad
  random_ben[[i]] <- preds_ben

  # Residuals
  samps <- sdmTMBextra::predict_mle_mcmc(m_sad, mcmc_iter = 201, mcmc_warmup = 200)
  mcmc_res <- residuals(m_sad, type = "mle-mcmc", mcmc_samples = samps)
  dd$res <- as.vector(mcmc_res)

  res_sad[[i]] <- dd

  samps <- sdmTMBextra::predict_mle_mcmc(m_ben, mcmc_iter = 201, mcmc_warmup = 200)
  mcmc_res <- residuals(m_ben, type = "mle-mcmc", mcmc_samples = samps)
  dd$res <- as.vector(mcmc_res)

  res_ben[[i]] <- dd


  # Ranges
  range_sad[[i]] <- tidy(m_sad, effects = "ran_pars") %>% filter(term == "range") %>% mutate(group = i, model = "saduria")
  range_ben[[i]] <- tidy(m_ben, effects = "ran_pars") %>% filter(term == "range") %>% mutate(group = i, model = "benthos")


  # Conditional effects: flounder
  nd_flounder <- data.frame(expand_grid(
    density_saduria_sc = c(quantile(d$density_saduria_sc, probs = 0.05),
                           quantile(d$density_saduria_sc, probs = 0.95)),
    flounder_density_sc = seq(quantile(dd$flounder_density_sc, probs = 0.05),
                              quantile(dd$flounder_density_sc, probs = 0.95),
                              length.out = 50))) %>%
    mutate(year = 2020,
           fyear = as.factor(2020),
           fquarter = as.factor(1),
           oxygen_sc = 0,
           depth_sc = 0,
           small_cod_density_sc = 0,
           fhaul_id = as.factor("2020_1_81")) # TODO: why do I need this when I specify re_form_iid!?

  preds_flounder_sad <- predict(m_sad, newdata = nd_flounder, re_form = NA, re_form_iid = NA, se_fit = TRUE)
  preds_flounder_ben <- predict(m_ben, newdata = nd_flounder, re_form = NA, re_form_iid = NA, se_fit = TRUE)

  pred_flounder_sad[[i]] <- preds_flounder_sad %>% mutate(group = i, xvar = "flounder")
  pred_flounder_ben[[i]] <- preds_flounder_ben %>% mutate(group = i, xvar = "flounder")

  # Conditional effects: cod
  nd_cod <- data.frame(expand_grid(
    density_saduria_sc = c(quantile(d$density_saduria_sc, probs = 0.05),
                           quantile(d$density_saduria_sc, probs = 0.95)),
    small_cod_density_sc = seq(quantile(dd$small_cod_density_sc, probs = 0.05),
                               quantile(dd$small_cod_density_sc, probs = 0.95),
                               length.out = 50))) %>%
    mutate(year = 2020,
           fyear = as.factor(2020),
           fquarter = as.factor(1),
           oxygen_sc = 0,
           depth_sc = 0,
           flounder_density_sc = 0,
           fhaul_id = as.factor("2020_1_81")) # TODO: why do I need this when I specify re_form_iid!?

  preds_cod_sad <- predict(m_sad, newdata = nd_cod, re_form = NA, re_form_iid = NA, se_fit = TRUE)
  preds_cod_ben <- predict(m_ben, newdata = nd_cod, re_form = NA, re_form_iid = NA, se_fit = TRUE)

  pred_cod_sad[[i]] <- preds_cod_sad %>% mutate(group = i, xvar = "cod")
  pred_cod_ben[[i]] <- preds_cod_ben %>% mutate(group = i, xvar = "cod")

  # Coefficients
  coefs_sad <- bind_rows(tidy(m_sad, effects = "fixed", conf.int = TRUE)) %>%
    mutate(species = "Cod (m)",
           response = "Saduria",
           sig = ifelse(estimate > 0 & conf.low > 0, "Y", "N"),
           sig = ifelse(estimate < 0 & conf.high < 0, "Y", sig))

  coefs_ben <- bind_rows(tidy(m_ben, effects = "fixed", conf.int = TRUE)) %>%
    mutate(species = "Cod (m)",
           response = "Saduria",
           sig = ifelse(estimate > 0 & conf.low > 0, "Y", "N"),
           sig = ifelse(estimate < 0 & conf.high < 0, "Y", sig))

  coef_sad[[i]] <- coefs_sad %>% mutate(group = i)
  coef_ben[[i]] <- coefs_ben %>% mutate(group = i)

}
filter: removed 5,765 rows (63%), 3,345 rows remaining
[1] "large cod"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
Spatiotemporal model fit by ML ['sdmTMB']
Formula: saduria_rel_weight ~ 0 + fyear + fquarter + depth_sc + oxygen_sc + 
 Formula:     small_cod_density_sc * density_saduria_sc + flounder_density_sc * 
 Formula:     density_saduria_sc
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: tweedie(link = 'log')
 
                                        coef.est coef.se
fyear2015                                  -8.36    1.02
fyear2016                                 -10.51    1.01
fyear2017                                 -10.77    0.95
fyear2018                                 -10.06    0.99
fyear2019                                 -11.32    1.14
fyear2020                                 -10.42    0.89
fyear2021                                 -10.29    0.90
fyear2022                                 -10.87    0.99
fquarter4                                  -0.46    0.29
depth_sc                                   -0.96    0.26
oxygen_sc                                  -0.02    0.22
small_cod_density_sc                        0.43    0.37
density_saduria_sc                          0.36    0.23
flounder_density_sc                        -0.36    0.25
small_cod_density_sc:density_saduria_sc     1.11    0.35
density_saduria_sc:flounder_density_sc     -0.15    0.28

Dispersion parameter: 0.14
Tweedie p: 1.47
Matérn range: 49.49
Spatial SD: 2.50
Spatiotemporal IID SD: 1.74
ML criterion at convergence: -721.647

See ?tidy.sdmTMB to extract these values as a data frame.
[1] "large cod"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
Spatiotemporal model fit by ML ['sdmTMB']
Formula: benthic_rel_weight ~ 0 + fyear + fquarter + depth_sc + small_cod_density_sc * 
 Formula:     oxygen_sc + flounder_density_sc * oxygen_sc
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: tweedie(link = 'log')
 
                               coef.est coef.se
fyear2015                         -6.56    0.30
fyear2016                         -6.56    0.26
fyear2017                         -6.76    0.22
fyear2018                         -6.37    0.24
fyear2019                         -7.25    0.29
fyear2020                         -6.75    0.22
fyear2021                         -6.44    0.21
fyear2022                         -6.59    0.24
fquarter4                          0.71    0.11
depth_sc                          -0.31    0.07
small_cod_density_sc               0.03    0.08
oxygen_sc                          0.16    0.06
flounder_density_sc               -0.22    0.06
small_cod_density_sc:oxygen_sc     0.10    0.09
oxygen_sc:flounder_density_sc      0.06    0.07

Dispersion parameter: 0.29
Tweedie p: 1.63
Matérn range: 28.79
Spatiotemporal IID SD: 0.92
ML criterion at convergence: -7117.386

See ?tidy.sdmTMB to extract these values as a data frame.
distinct: removed 3,057 rows (91%), 288 rows remaining

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.002759 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 27.59 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 17.021 seconds (Warm-up)
Chain 1:                0.051 seconds (Sampling)
Chain 1:                17.072 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.002573 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 25.73 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 10.564 seconds (Warm-up)
Chain 1:                0.024 seconds (Sampling)
Chain 1:                10.588 seconds (Total)
Chain 1: 
filter: removed 4 rows (80%), one row remaining
mutate: new variable 'group' (character) with one unique value and 0% NA
        new variable 'model' (character) with one unique value and 0% NA
filter: removed 3 rows (75%), one row remaining
mutate: new variable 'group' (character) with one unique value and 0% NA
        new variable 'model' (character) with one unique value and 0% NA
mutate: new variable 'year' (double) with one unique value and 0% NA
        new variable 'fyear' (factor) with one unique value and 0% NA
        new variable 'fquarter' (factor) with one unique value and 0% NA
        new variable 'oxygen_sc' (double) with one unique value and 0% NA
        new variable 'depth_sc' (double) with one unique value and 0% NA
        new variable 'small_cod_density_sc' (double) with one unique value and 0% NA
        new variable 'fhaul_id' (factor) with one unique value and 0% NA
mutate: new variable 'group' (character) with one unique value and 0% NA
        new variable 'xvar' (character) with one unique value and 0% NA
mutate: new variable 'group' (character) with one unique value and 0% NA
        new variable 'xvar' (character) with one unique value and 0% NA
mutate: new variable 'year' (double) with one unique value and 0% NA
        new variable 'fyear' (factor) with one unique value and 0% NA
        new variable 'fquarter' (factor) with one unique value and 0% NA
        new variable 'oxygen_sc' (double) with one unique value and 0% NA
        new variable 'depth_sc' (double) with one unique value and 0% NA
        new variable 'flounder_density_sc' (double) with one unique value and 0% NA
        new variable 'fhaul_id' (factor) with one unique value and 0% NA
mutate: new variable 'group' (character) with one unique value and 0% NA
        new variable 'xvar' (character) with one unique value and 0% NA
mutate: new variable 'group' (character) with one unique value and 0% NA
        new variable 'xvar' (character) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
        new variable 'response' (character) with one unique value and 0% NA
        new variable 'sig' (character) with 2 unique values and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
        new variable 'response' (character) with one unique value and 0% NA
        new variable 'sig' (character) with 2 unique values and 0% NA
mutate: new variable 'group' (character) with one unique value and 0% NA
mutate: new variable 'group' (character) with one unique value and 0% NA
filter: removed 7,196 rows (79%), 1,914 rows remaining
[1] "small cod"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
Spatiotemporal model fit by ML ['sdmTMB']
Formula: saduria_rel_weight ~ 0 + fyear + fquarter + depth_sc + oxygen_sc + 
 Formula:     small_cod_density_sc * density_saduria_sc + flounder_density_sc * 
 Formula:     density_saduria_sc
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: tweedie(link = 'log')
 
                                        coef.est coef.se
fyear2015                                  -8.16    1.20
fyear2016                                  -9.30    1.14
fyear2017                                  -9.99    1.11
fyear2018                                  -9.72    1.28
fyear2019                                  -8.40    1.24
fyear2020                                  -9.12    1.09
fyear2021                                  -9.75    1.08
fyear2022                                  -9.14    1.10
fquarter4                                  -1.34    0.34
depth_sc                                   -0.73    0.27
oxygen_sc                                  -0.47    0.23
small_cod_density_sc                        0.38    0.34
density_saduria_sc                          0.55    0.24
flounder_density_sc                        -0.77    0.34
small_cod_density_sc:density_saduria_sc     0.25    0.30
density_saduria_sc:flounder_density_sc     -0.15    0.39

Dispersion parameter: 0.19
Tweedie p: 1.50
Matérn range: 87.59
Spatial SD: 2.26
Spatiotemporal IID SD: 1.33
ML criterion at convergence: -591.883

See ?tidy.sdmTMB to extract these values as a data frame.
[1] "small cod"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
Spatiotemporal model fit by ML ['sdmTMB']
Formula: benthic_rel_weight ~ 0 + fyear + fquarter + depth_sc + small_cod_density_sc * 
 Formula:     oxygen_sc + flounder_density_sc * oxygen_sc
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: tweedie(link = 'log')
 
                               coef.est coef.se
fyear2015                         -5.88    0.29
fyear2016                         -5.85    0.25
fyear2017                         -5.74    0.21
fyear2018                         -5.65    0.24
fyear2019                         -6.22    0.32
fyear2020                         -5.73    0.21
fyear2021                         -5.94    0.19
fyear2022                         -5.66    0.21
fquarter4                          0.42    0.10
depth_sc                          -0.34    0.07
small_cod_density_sc               0.07    0.07
oxygen_sc                          0.16    0.06
flounder_density_sc               -0.04    0.05
small_cod_density_sc:oxygen_sc    -0.20    0.10
oxygen_sc:flounder_density_sc      0.02    0.07

Dispersion parameter: 0.11
Tweedie p: 1.54
Matérn range: 33.19
Spatiotemporal IID SD: 0.72
ML criterion at convergence: -5129.440

See ?tidy.sdmTMB to extract these values as a data frame.
distinct: removed 1,647 rows (86%), 267 rows remaining

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.002716 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 27.16 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 41.203 seconds (Warm-up)
Chain 1:                0.393 seconds (Sampling)
Chain 1:                41.596 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.002803 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 28.03 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 9.571 seconds (Warm-up)
Chain 1:                0.022 seconds (Sampling)
Chain 1:                9.593 seconds (Total)
Chain 1: 
filter: removed 4 rows (80%), one row remaining
mutate: new variable 'group' (character) with one unique value and 0% NA
        new variable 'model' (character) with one unique value and 0% NA
filter: removed 3 rows (75%), one row remaining
mutate: new variable 'group' (character) with one unique value and 0% NA
        new variable 'model' (character) with one unique value and 0% NA
mutate: new variable 'year' (double) with one unique value and 0% NA
        new variable 'fyear' (factor) with one unique value and 0% NA
        new variable 'fquarter' (factor) with one unique value and 0% NA
        new variable 'oxygen_sc' (double) with one unique value and 0% NA
        new variable 'depth_sc' (double) with one unique value and 0% NA
        new variable 'small_cod_density_sc' (double) with one unique value and 0% NA
        new variable 'fhaul_id' (factor) with one unique value and 0% NA
mutate: new variable 'group' (character) with one unique value and 0% NA
        new variable 'xvar' (character) with one unique value and 0% NA
mutate: new variable 'group' (character) with one unique value and 0% NA
        new variable 'xvar' (character) with one unique value and 0% NA
mutate: new variable 'year' (double) with one unique value and 0% NA
        new variable 'fyear' (factor) with one unique value and 0% NA
        new variable 'fquarter' (factor) with one unique value and 0% NA
        new variable 'oxygen_sc' (double) with one unique value and 0% NA
        new variable 'depth_sc' (double) with one unique value and 0% NA
        new variable 'flounder_density_sc' (double) with one unique value and 0% NA
        new variable 'fhaul_id' (factor) with one unique value and 0% NA
mutate: new variable 'group' (character) with one unique value and 0% NA
        new variable 'xvar' (character) with one unique value and 0% NA
mutate: new variable 'group' (character) with one unique value and 0% NA
        new variable 'xvar' (character) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
        new variable 'response' (character) with one unique value and 0% NA
        new variable 'sig' (character) with 2 unique values and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
        new variable 'response' (character) with one unique value and 0% NA
        new variable 'sig' (character) with 2 unique values and 0% NA
mutate: new variable 'group' (character) with one unique value and 0% NA
mutate: new variable 'group' (character) with one unique value and 0% NA
filter: removed 5,259 rows (58%), 3,851 rows remaining
[1] "flounder"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
Spatiotemporal model fit by ML ['sdmTMB']
Formula: saduria_rel_weight ~ 0 + fyear + fquarter + depth_sc + oxygen_sc + 
 Formula:     small_cod_density_sc * density_saduria_sc + flounder_density_sc * 
 Formula:     density_saduria_sc
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: tweedie(link = 'log')
 
                                        coef.est coef.se
fyear2015                                  -8.33    1.60
fyear2016                                  -6.61    1.43
fyear2017                                  -9.24    1.42
fyear2018                                  -8.66    1.49
fyear2019                                  -8.84    1.53
fyear2020                                  -8.90    1.44
fyear2021                                 -13.05    1.79
fyear2022                                  -9.64    1.47
fquarter4                                  -0.46    0.20
depth_sc                                   -0.43    0.18
oxygen_sc                                  -0.02    0.21
small_cod_density_sc                        0.10    0.26
density_saduria_sc                         -0.03    0.17
flounder_density_sc                        -0.80    0.23
small_cod_density_sc:density_saduria_sc    -0.01    0.30
density_saduria_sc:flounder_density_sc     -0.68    0.21

Dispersion parameter: 0.19
Tweedie p: 1.49
Matérn range: 111.89
Spatiotemporal IID SD: 3.08
ML criterion at convergence: -1764.896

See ?tidy.sdmTMB to extract these values as a data frame.
[1] "flounder"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
Spatiotemporal model fit by ML ['sdmTMB']
Formula: benthic_rel_weight ~ 0 + fyear + fquarter + depth_sc + small_cod_density_sc * 
 Formula:     oxygen_sc + flounder_density_sc * oxygen_sc
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: tweedie(link = 'log')
 
                               coef.est coef.se
fyear2015                         -5.01    0.34
fyear2016                         -5.17    0.28
fyear2017                         -5.15    0.26
fyear2018                         -5.59    0.30
fyear2019                         -5.27    0.30
fyear2020                         -4.68    0.25
fyear2021                         -5.40    0.29
fyear2022                         -5.02    0.27
fquarter4                          0.11    0.09
depth_sc                          -0.22    0.09
small_cod_density_sc              -0.03    0.06
oxygen_sc                         -0.02    0.06
flounder_density_sc               -0.11    0.05
small_cod_density_sc:oxygen_sc    -0.15    0.09
oxygen_sc:flounder_density_sc     -0.01    0.06

Dispersion parameter: 0.18
Tweedie p: 1.51
Matérn range: 29.35
Spatial SD: 1.20
Spatiotemporal IID SD: 0.75
ML criterion at convergence: -5507.909

See ?tidy.sdmTMB to extract these values as a data frame.
distinct: removed 3,590 rows (93%), 261 rows remaining

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.002122 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 21.22 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 34.591 seconds (Warm-up)
Chain 1:                0.157 seconds (Sampling)
Chain 1:                34.748 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.002502 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 25.02 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 14.076 seconds (Warm-up)
Chain 1:                0.047 seconds (Sampling)
Chain 1:                14.123 seconds (Total)
Chain 1: 
filter: removed 3 rows (75%), one row remaining
mutate: new variable 'group' (character) with one unique value and 0% NA
        new variable 'model' (character) with one unique value and 0% NA
filter: removed 4 rows (80%), one row remaining
mutate: new variable 'group' (character) with one unique value and 0% NA
        new variable 'model' (character) with one unique value and 0% NA
mutate: new variable 'year' (double) with one unique value and 0% NA
        new variable 'fyear' (factor) with one unique value and 0% NA
        new variable 'fquarter' (factor) with one unique value and 0% NA
        new variable 'oxygen_sc' (double) with one unique value and 0% NA
        new variable 'depth_sc' (double) with one unique value and 0% NA
        new variable 'small_cod_density_sc' (double) with one unique value and 0% NA
        new variable 'fhaul_id' (factor) with one unique value and 0% NA
mutate: new variable 'group' (character) with one unique value and 0% NA
        new variable 'xvar' (character) with one unique value and 0% NA
mutate: new variable 'group' (character) with one unique value and 0% NA
        new variable 'xvar' (character) with one unique value and 0% NA
mutate: new variable 'year' (double) with one unique value and 0% NA
        new variable 'fyear' (factor) with one unique value and 0% NA
        new variable 'fquarter' (factor) with one unique value and 0% NA
        new variable 'oxygen_sc' (double) with one unique value and 0% NA
        new variable 'depth_sc' (double) with one unique value and 0% NA
        new variable 'flounder_density_sc' (double) with one unique value and 0% NA
        new variable 'fhaul_id' (factor) with one unique value and 0% NA
mutate: new variable 'group' (character) with one unique value and 0% NA
        new variable 'xvar' (character) with one unique value and 0% NA
mutate: new variable 'group' (character) with one unique value and 0% NA
        new variable 'xvar' (character) with one unique value and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
        new variable 'response' (character) with one unique value and 0% NA
        new variable 'sig' (character) with 2 unique values and 0% NA
mutate: new variable 'species' (character) with one unique value and 0% NA
        new variable 'response' (character) with one unique value and 0% NA
        new variable 'sig' (character) with 2 unique values and 0% NA
mutate: new variable 'group' (character) with one unique value and 0% NA
mutate: new variable 'group' (character) with one unique value and 0% NA

Now do a separate model for adult cod looking at total prey

  dd <- filter(d, group == "large cod")
filter: removed 5,765 rows (63%), 3,345 rows remaining
  mesh <- make_mesh(dd,
                    xy_cols = c("X", "Y"),
                    cutoff = 10)

  # Total model
  # NOTE: turning off spatial here due to convergence and AIC
  m_tot <- sdmTMB(tot_rel_weight ~ 0 + fyear + fquarter + depth_sc + #(1|fhaul_id) +
                    large_cod_density_sc*oxygen_sc + flounder_density_sc*oxygen_sc,
                  data = dd,
                  mesh = mesh,
                  family = tweedie(),
                  spatiotemporal = "IID", 
                  spatial = "off",
                  time = "year")
  
  sanity(m_tot)
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
  print(m_tot)
Spatiotemporal model fit by ML ['sdmTMB']
Formula: tot_rel_weight ~ 0 + fyear + fquarter + depth_sc + large_cod_density_sc * 
 Formula:     oxygen_sc + flounder_density_sc * oxygen_sc
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: tweedie(link = 'log')
 
                               coef.est coef.se
fyear2015                         -4.20    0.20
fyear2016                         -4.69    0.17
fyear2017                         -4.60    0.14
fyear2018                         -4.32    0.14
fyear2019                         -4.75    0.21
fyear2020                         -4.57    0.13
fyear2021                         -4.39    0.13
fyear2022                         -4.67    0.15
fquarter4                         -0.19    0.09
depth_sc                          -0.04    0.05
large_cod_density_sc               0.13    0.06
oxygen_sc                         -0.03    0.05
flounder_density_sc               -0.03    0.04
large_cod_density_sc:oxygen_sc     0.10    0.11
oxygen_sc:flounder_density_sc     -0.07    0.05

Dispersion parameter: 0.62
Tweedie p: 1.71
Matérn range: 15.85
Spatiotemporal IID SD: 0.73
ML criterion at convergence: -7227.381

See ?tidy.sdmTMB to extract these values as a data frame.
  # Residuals
  samps <- sdmTMBextra::predict_mle_mcmc(m_tot, mcmc_iter = 201, mcmc_warmup = 200)

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.002636 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 26.36 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 7.555 seconds (Warm-up)
Chain 1:                0.024 seconds (Sampling)
Chain 1:                7.579 seconds (Total)
Chain 1: 
  mcmc_res <- residuals(m_tot, type = "mle-mcmc", mcmc_samples = samps)
  dd$res <- as.vector(mcmc_res)
    
  res_tot <- dd
  
  # Range
  range_tot <- tidy(m_tot, effects = "ran_pars") %>% filter(term == "range") %>% mutate(group = "large cod", model = "total")
filter: removed 3 rows (75%), one row remaining
mutate: new variable 'group' (character) with one unique value and 0% NA
        new variable 'model' (character) with one unique value and 0% NA
  # Spatial and spatiotemporal random effects
  d_haul <- dd %>% 
    distinct(haul_id, .keep_all = TRUE)
distinct: removed 3,057 rows (91%), 288 rows remaining
  preds_tot <- predict(m_tot, newdata = d_haul)
  
  # Coefficients
  coefs_tot <- bind_rows(tidy(m_tot, effects = "fixed", conf.int = TRUE)) %>% 
    mutate(species = "Cod (m)",
           response = "Saduria",
           sig = ifelse(estimate > 0 & conf.low > 0, "Y", "N"),
           sig = ifelse(estimate < 0 & conf.high < 0, "Y", sig))
mutate: new variable 'species' (character) with one unique value and 0% NA
        new variable 'response' (character) with one unique value and 0% NA
        new variable 'sig' (character) with 2 unique values and 0% NA
  coefs_tot <- coefs_tot %>% mutate(group = "large cod")
mutate: new variable 'group' (character) with one unique value and 0% NA

Make dataframes

coef_df <- bind_rows(bind_rows(coef_sad) %>% mutate(model = "Saduria"),
                     bind_rows(coef_ben) %>% mutate(model = "Benthos"))
mutate: new variable 'model' (character) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NA
coef_df <- coef_df %>% bind_rows(coefs_tot %>% mutate(model = "Total"))
mutate: new variable 'model' (character) with one unique value and 0% NA
pred_cod_df <- bind_rows(bind_rows(pred_cod_sad) %>% mutate(model = "Saduria"),
                         bind_rows(pred_cod_ben) %>% mutate(model = "Benthos"))
mutate: new variable 'model' (character) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NA
pred_flounder_df <- bind_rows(bind_rows(pred_flounder_sad) %>% mutate(model = "Saduria"),
                              bind_rows(pred_flounder_ben) %>% mutate(model = "Benthos"))
mutate: new variable 'model' (character) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NA
res_df <- bind_rows(bind_rows(res_sad) %>% mutate(model = "Saduria"),
                    bind_rows(res_ben) %>% mutate(model = "Benthos"))
mutate: new variable 'model' (character) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NA
res_df <- res_df %>% bind_rows(res_tot %>% mutate(model = "Total"))
mutate: new variable 'model' (character) with one unique value and 0% NA
random_df <- bind_rows(bind_rows(random_sad) %>% mutate(model = "Saduria"),
                       bind_rows(random_ben) %>% mutate(model = "Benthos"))
mutate: new variable 'model' (character) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NA
random_df <- random_df %>% bind_rows(preds_tot %>% mutate(model = "Total"))
mutate: new variable 'model' (character) with one unique value and 0% NA
range_df <- bind_rows(range_tot, bind_rows(range_ben), bind_rows(range_sad))

Plot spatial random effects

random_df <- random_df %>%
  mutate(group = str_to_sentence(group))
mutate: changed 1,920 values (100%) of 'group' (0 new NA)
# Saduria
sad_omega <- plot_map_fc +
  geom_point(data = random_df %>% filter(model == "Saduria" & !group == "Flounder"), aes(X*1000, Y*1000, color = omega_s), size = 0.9) +
  scale_color_gradient2() +
  facet_wrap(~ group) +
  labs(color = "Spatial\nrandom effect") +
  theme(axis.text.x = element_text(angle = 90),
        legend.position = "right",
        legend.direction = "vertical",
        legend.key.width = unit(0.4, "cm"),
        legend.key.height = unit(0.4, "cm"))
filter: removed 1,365 rows (71%), 555 rows remaining
ggsave(paste0(home, "/figures/supp/omega_sad.pdf"), width = 17, height = 7, units = "cm")


# Now do benthos (only for flounder)
ben_omega <- plot_map_fc +
  geom_point(data = random_df %>% filter(model == "Benthos" & group == "Flounder"), aes(X*1000, Y*1000, color = omega_s), size = 0.9) +
  scale_color_gradient2() +
  facet_wrap(~ group) +
  labs(color = "Spatial\nrandom effect") +
  theme(axis.text.x = element_text(angle = 90),
        legend.position = "right",
        legend.direction = "vertical",
        legend.key.width = unit(0.4, "cm"),
        legend.key.height = unit(0.4, "cm"))
filter: removed 1,659 rows (86%), 261 rows remaining
ben_omega

ggsave(paste0(home, "/figures/supp/omega_ben.pdf"), width = 11, height = 11, units = "cm")

Plot spatiotemporal random effects

# Saduria
sad_eps_sc <- plot_map_fc +
  geom_point(data = random_df %>% filter(model == "Saduria" & group == "Small cod"), aes(X*1000, Y*1000, color = epsilon_st), size = 0.9) +
  scale_color_gradient2() +
  facet_wrap(~ year) +
  labs(color = "Spatiotemporal\nrandom effect") +
  theme(legend.position = c(0.84, 0.16),
        axis.text.x = element_text(angle = 90))
filter: removed 1,653 rows (86%), 267 rows remaining
sad_eps_sc

ggsave(paste0(home, "/figures/supp/epsilon_sad_small_cod.pdf"), width = 17, height = 17, units = "cm")

sad_eps_lc <- plot_map_fc +
  geom_point(data = random_df %>% filter(model == "Saduria" & group == "Large cod"), aes(X*1000, Y*1000, color = epsilon_st), size = 0.9) +
  scale_color_gradient2() +
  facet_wrap(~ year) +
  labs(color = "Spatiotemporal\nrandom effect") +
  theme(legend.position = c(0.84, 0.16),
        axis.text.x = element_text(angle = 90))
filter: removed 1,632 rows (85%), 288 rows remaining
sad_eps_lc

ggsave(paste0(home, "/figures/supp/epsilon_sad_large_cod.pdf"), width = 17, height = 17, units = "cm")

sad_eps_f <- plot_map_fc +
  geom_point(data = random_df %>% filter(model == "Saduria" & group == "Flounder"), aes(X*1000, Y*1000, color = epsilon_st), size = 0.9) +
  scale_color_gradient2() +
  facet_wrap(~ year) +
  labs(color = "Spatiotemporal\nrandom effect") +
  theme(legend.position = c(0.84, 0.16),
        axis.text.x = element_text(angle = 90))
filter: removed 1,659 rows (86%), 261 rows remaining
sad_eps_f

ggsave(paste0(home, "/figures/supp/epsilon_sad_flounder.pdf"), width = 17, height = 17, units = "cm")


# Benthos
ben_eps_sc <- plot_map_fc +
  geom_point(data = random_df %>% filter(model == "Benthos" & group == "Small cod"), aes(X*1000, Y*1000, color = epsilon_st), size = 0.9) +
  scale_color_gradient2() +
  facet_wrap(~ year) +
  labs(color = "Spatiotemporal\nrandom effect") +
  theme(legend.position = c(0.84, 0.16),
        axis.text.x = element_text(angle = 90))
filter: removed 1,653 rows (86%), 267 rows remaining
ben_eps_sc

ggsave(paste0(home, "/figures/supp/epsilon_ben_small_cod.pdf"), width = 17, height = 17, units = "cm")

ben_eps_lc <- plot_map_fc +
  geom_point(data = random_df %>% filter(model == "Benthos" & group == "Large cod"), aes(X*1000, Y*1000, color = epsilon_st), size = 0.9) +
  scale_color_gradient2() +
  facet_wrap(~ year) +
  labs(color = "Spatiotemporal\nrandom effect") +
  theme(legend.position = c(0.84, 0.16),
        axis.text.x = element_text(angle = 90))
filter: removed 1,632 rows (85%), 288 rows remaining
ben_eps_lc

ggsave(paste0(home, "/figures/supp/epsilon_ben_large_cod.pdf"), width = 17, height = 17, units = "cm")

ben_eps_f <- plot_map_fc +
  geom_point(data = random_df %>% filter(model == "Benthos" & group == "Flounder"), aes(X*1000, Y*1000, color = epsilon_st), size = 0.9) +
  scale_color_gradient2() +
  facet_wrap(~ year) +
  labs(color = "Spatiotemporal\nrandom effect") +
  theme(legend.position = c(0.84, 0.16),
        axis.text.x = element_text(angle = 90))
filter: removed 1,659 rows (86%), 261 rows remaining
ben_eps_f

ggsave(paste0(home, "/figures/supp/epsilon_ben_flounder.pdf"), width = 17, height = 17, units = "cm")


# Total
tot_eps <- plot_map_fc + 
    geom_point(data = preds_tot, aes(X*1000, Y*1000, color = epsilon_st), size = 0.9) + 
    scale_color_gradient2() +
    facet_wrap(~ year) +
    labs(color = "Spatiotemporal\nrandom effect") +
    theme(legend.position = c(0.84, 0.16),
          axis.text.x = element_text(angle = 90))

tot_eps

Plot range

pal <- brewer.pal(n = 8, name = "Dark2")[c(2, 7, 6)]

range_df %>% 
  mutate(group = str_to_sentence(group),
         model = str_to_sentence(model)) %>% 
  ggplot(aes(model, estimate, color = group)) + 
  geom_point(size = 2) +
  geom_hline(yintercept = 15, linetype = 2, alpha = 0.5) +
  scale_color_manual(values = pal) + 
  labs(x = "Prey group", y = "Range (km)", color = "Group") + 
  theme(aspect.ratio = 1,
        legend.position = c(0.86, 0.86)) 
mutate: changed 7 values (100%) of 'group' (0 new NA)
        changed 7 values (100%) of 'model' (0 new NA)

ggsave(paste0(home, "/figures/supp/ranges.pdf"), width = 11, height = 11, units = "cm")

Plot residuals

# Plot residuals
res_df |> 
  mutate(group = str_to_title(group)) |> 
  ggplot(aes(sample = res)) +
  stat_qq(size = 0.75, shape = 21, fill = NA) +
  facet_grid(model ~ group) +
  stat_qq_line() +
  labs(y = "Sample Quantiles", x = "Theoretical Quantiles") +
  theme(aspect.ratio = 1)
mutate: changed 21,565 values (100%) of 'group' (0 new NA)

ggsave(paste0(home, "/figures/supp/qq_relative_prey_weight.pdf"), width = 17, height = 17, units = "cm")

Plot coefficients

coef_df$term %>% unique()
 [1] "fyear2015"                              
 [2] "fyear2016"                              
 [3] "fyear2017"                              
 [4] "fyear2018"                              
 [5] "fyear2019"                              
 [6] "fyear2020"                              
 [7] "fyear2021"                              
 [8] "fyear2022"                              
 [9] "fquarter4"                              
[10] "depth_sc"                               
[11] "oxygen_sc"                              
[12] "small_cod_density_sc"                   
[13] "density_saduria_sc"                     
[14] "flounder_density_sc"                    
[15] "small_cod_density_sc:density_saduria_sc"
[16] "density_saduria_sc:flounder_density_sc" 
[17] "small_cod_density_sc:oxygen_sc"         
[18] "oxygen_sc:flounder_density_sc"          
[19] "large_cod_density_sc"                   
[20] "large_cod_density_sc:oxygen_sc"         
pal <- brewer.pal(n = 8, name = "Dark2")[c(2, 7, 6)]

# Fix some names
coef_df2 <- coef_df %>%
  filter(!grepl('year', term)) %>%
  filter(!grepl('quarter', term)) %>%
  mutate(term = str_remove_all(term, "_sc"),
         term = str_remove_all(term, "density"),
         term = str_replace_all(term, "_", ""),
         term = str_replace_all(term, "geco", "ge co"),
         term = str_replace_all(term, "llco", "ll co"),
         term = str_replace(term, ":", " × "),
         term = str_to_sentence(term),
         group = str_to_sentence(group))
filter: removed 56 rows (52%), 52 rows remaining
filter: removed 7 rows (13%), 45 rows remaining
mutate: changed 45 values (100%) of 'term' (0 new NA)
        changed 45 values (100%) of 'group' (0 new NA)
ggplot(coef_df2, aes(estimate, term, color = group, alpha = sig)) +
  geom_stripped_rows(aes(y = term), inherit.aes = FALSE) +
  facet_wrap2(~model, ncol = 2, scales = "free") +
  geom_vline(xintercept = 0, linetype = 2, alpha = 0.5, color = "gray10", linewidth = 0.2) +
  geom_point(position = position_dodge(width = 0.5), size = 1.5) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), width = 0,
                position = position_dodge(width = 0.5)) +
  scale_alpha_manual(values = c(0.4, 1)) +
  scale_color_manual(values = pal) +
  labs(x = "", y = "Standardized coefficient", alpha = "95% CI crossing 0", color = "Group") +
  guides(color = guide_legend(title.position = "top", title.hjust = 0.5),
         alpha = guide_legend(title.position = "top", title.hjust = 0.5)) +
  theme(legend.position = c(0.75, 0.2),
        legend.direction = "vertical",
        legend.box.spacing = unit(-3, "pt"),
        legend.margin = margin(0, 0, 0, 0)) +
  NULL

ggsave(paste0(home, "/figures/coefs.pdf"), width = 17, height = 13, units = "cm")

Plot year and quarter coefficients

# Fix some names
coef_df3 <- coef_df %>%
  filter(grepl('year', term)) %>%
  mutate(term = str_remove_all(term, "fyear"),
         group = str_to_sentence(group),
         term = as.numeric(term))
filter: removed 52 rows (48%), 56 rows remaining
mutate: converted 'term' from character to double (0 new NA)
        changed 56 values (100%) of 'group' (0 new NA)
ggplot(coef_df3, aes(term, estimate, color = group, fill = group)) +
  facet_wrap(~model, scales = "free", ncol = 1) +
  geom_line(position = position_dodge(width = 0.5)) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.2, color = NA) +
  scale_color_manual(values = pal) +
  scale_fill_manual(values = pal) +
  labs(x = "Year", y = "Standardized coefficient", color = "") +
  guides(color = guide_legend(title.position = "top", title.hjust = 0.5, ncol = 3),
         fill = "none") +
  theme(legend.position = c(0.5, 0.99),
        legend.direction = "vertical",
        legend.box.spacing = unit(-3, "pt"),
        legend.margin = margin(0, 0, 0, 0),
        strip.text.x.top = element_text(angle = 0, hjust = 0)) +
  NULL

ggsave(paste0(home, "/figures/supp/coefs_year.pdf"), width = 11, height = 21, units = "cm")


# Now do quarter
coef_df5 <- coef_df %>%
  filter(term %in% c("fquarter4")) %>% 
  mutate(group = str_to_sentence(group))
filter: removed 101 rows (94%), 7 rows remaining
mutate: changed 7 values (100%) of 'group' (0 new NA)
ggplot(coef_df5, aes(estimate, model, color = group, alpha = sig)) +
  geom_vline(xintercept = 0, linetype = 2, alpha = 0.5, color = "gray10", linewidth = 0.2) +
  geom_point(position = position_dodge(width = 0.5), size = 1.5) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), width = 0,
                position = position_dodge(width = 0.5)) +
  scale_alpha_manual(values = c(0.4, 1)) +
  scale_color_manual(values = pal) +
  labs(x = "", y = "Quarter 4 effect", alpha = "95% CI crossing 0", color = "Group") +
  guides(color = guide_legend(title.position = "top", title.hjust = 0.5),
         alpha = guide_legend(title.position = "top", title.hjust = 0.5)) +
  theme(legend.position = "bottom",
        legend.direction = "horizontal",
        legend.box = "horizontal",
        legend.box.spacing = unit(-3, "pt"),
        legend.margin = margin(0, 0, 0, 0))

ggsave(paste0(home, "/figures/supp/coefs_quarter.pdf"), width = 17, height = 11, units = "cm")

Conditional effects

# Which CI?
# https://www.calculator.net/confidence-interval-calculator.html
pred_df <- bind_rows(pred_cod_df, pred_flounder_df) %>%
  mutate(group = str_to_sentence(group),
         sad = ifelse(density_saduria_sc == min(density_saduria_sc), "Low", "High"))
mutate: changed 1,200 values (100%) of 'group' (0 new NA)
        new variable 'sad' (character) with 2 unique values and 0% NA
# 75% CI!!

# x small cod for different saduria (y = small cod)
# x flounder cod for different saduria (y = all)

ggplot(pred_df %>% filter(model == "Saduria" & xvar == "cod" & group == "Small cod"),
       aes(small_cod_density_sc, exp(est), color = sad, fill = sad)) +
  geom_ribbon(aes(ymin = exp(est - 1.036*est_se), ymax = exp(est + 1.036*est_se)),
             alpha = 0.3, color = NA) +
  geom_line() +
  #coord_cartesian(ylim = c(0, 0.0001)) +
  facet_wrap(~group, scales = "free") +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  labs(x = "Small cod density", y = "Relative saduria weight",
       color = "Saduria", fill = "Saduria") +
  theme(legend.position = c(0.88, 0.88))
filter: removed 1,100 rows (92%), 100 rows remaining

#ggsave(paste0(home, "/figures/conditional_saduria_scod.pdf"), width = 11, height = 11, units = "cm")

ggplot(pred_df %>% filter(model == "Saduria" & xvar == "flounder" & !group == "Large cod"),
       aes(flounder_density_sc, exp(est), color = sad, fill = sad)) +
  geom_ribbon(aes(ymin = exp(est - 1.150*est_se), ymax = exp(est + 1.150*est_se)),
              alpha = 0.3, color = NA) +
  geom_line() +
  facet_wrap(~group, scales = "free", ncol = 2) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  coord_cartesian(xlim = c(-0.5, 0.5)) +
  labs(x = "Flounder density", y = "Relative saduria weight",
       color = "Saduria", fill = "Saduria") + 
  theme(legend.position = c(0.95, 0.88),
        strip.text.x.top = element_text(angle = 0, hjust = 0))
filter: removed 1,000 rows (83%), 200 rows remaining

ggsave(paste0(home, "/figures/conditional_saduria_flounder.pdf"), width = 17, height = 9, units = "cm")

# Conditional effects of oxygen and flounder for the benthos model

  nd_flounder_oxy <- data.frame(expand_grid(
    oxygen_sc = c(quantile(d$oxygen_sc, probs = 0.05),
                  quantile(d$oxygen_sc, probs = 0.95)),
    flounder_density_sc = seq(quantile(dd$flounder_density_sc, probs = 0.05),
                              quantile(dd$flounder_density_sc, probs = 0.95),
                              length.out = 50))) %>%
    mutate(year = 2020,
           fyear = as.factor(2020),
           fquarter = as.factor(1),
           density_saduria_sc = 0, 
           depth_sc = 0,
           small_cod_density_sc = 0,
           fhaul_id = as.factor("2020_1_81")) # TODO: why do I need this when I specify re_form_iid!?
mutate: new variable 'year' (double) with one unique value and 0% NA
        new variable 'fyear' (factor) with one unique value and 0% NA
        new variable 'fquarter' (factor) with one unique value and 0% NA
        new variable 'density_saduria_sc' (double) with one unique value and 0% NA
        new variable 'depth_sc' (double) with one unique value and 0% NA
        new variable 'small_cod_density_sc' (double) with one unique value and 0% NA
        new variable 'fhaul_id' (factor) with one unique value and 0% NA
  preds_flounder_oxy_ben <- predict(m_ben, newdata = nd_flounder_oxy, re_form = NA, re_form_iid = NA, se_fit = TRUE)

  
  ggplot(preds_flounder_oxy_ben, aes(flounder_density_sc, exp(est), color = factor(round(oxygen_sc)), fill = factor(round(oxygen_sc)))) +
    geom_ribbon(aes(ymin = exp(est - 1.150*est_se), ymax = exp(est + 1.150*est_se)),
                alpha = 0.3, color = NA) +
    geom_line() +
    #facet_wrap(~group, scales = "free", ncol = 2) +
    scale_color_brewer(palette = "Dark2") +
    scale_fill_brewer(palette = "Dark2") +
    coord_cartesian(xlim = c(-0.5, 0.5)) +
    labs(x = "Flounder density", y = "Relative saduria weight",
         color = "Oxygen", fill = "Oxygen") + 
    theme(legend.position = c(0.95, 0.88),
          strip.text.x.top = element_text(angle = 0, hjust = 0))

Showing conditional effects of oxygen on small cod feeding on benthos

  dd <- filter(d, group == "small cod")
filter: removed 7,196 rows (79%), 1,914 rows remaining
  mesh <- make_mesh(dd,
                    xy_cols = c("X", "Y"),
                    cutoff = 10)

  # Benthic model
  m_ben <- sdmTMB(benthic_rel_weight ~ 0 + fyear + fquarter + depth_sc + #(1|fhaul_id) +
                    small_cod_density_sc*oxygen_sc + flounder_density_sc*oxygen_sc,
                  data = dd,
                  mesh = mesh,
                  family = tweedie(),
                  spatiotemporal = "IID", 
                  spatial = "off",
                  time = "year")

  sanity(m_ben)
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
  nd <- data.frame(oxygen = seq(quantile(d$oxygen, probs = 0.05), quantile(d$oxygen, probs = 0.95),
                                length.out = 50)) %>%
    mutate(year = 2020,
           fyear = as.factor(2020),
           fquarter = as.factor(1),
           density_saduria_sc = 0,
           flounder_density_sc = 0, 
           depth_sc = 0,
           small_cod_density_sc = 0,
           fhaul_id = as.factor("2020_1_81")) %>% 
    mutate(oxygen_sc = (oxygen - mean(d$oxygen)) / sd(d$oxygen)) 
mutate: new variable 'year' (double) with one unique value and 0% NA
        new variable 'fyear' (factor) with one unique value and 0% NA
        new variable 'fquarter' (factor) with one unique value and 0% NA
        new variable 'density_saduria_sc' (double) with one unique value and 0% NA
        new variable 'flounder_density_sc' (double) with one unique value and 0% NA
        new variable 'depth_sc' (double) with one unique value and 0% NA
        new variable 'small_cod_density_sc' (double) with one unique value and 0% NA
        new variable 'fhaul_id' (factor) with one unique value and 0% NA
mutate: new variable 'oxygen_sc' (double) with 50 unique values and 0% NA
  p <- predict(m_ben, newdata = nd, re_form = NA, re_form_iid = NA, se_fit = TRUE)
  
  ggplot(p, aes(oxygen, exp(est))) + 
    geom_line() + 
    theme_sleek(base_size = 14) + 
    geom_hline(yintercept = 0.0025, col = "red") +
    geom_hline(yintercept = 0.004, col = "red") +
    geom_vline(xintercept = 4, col = "red") +
    geom_vline(xintercept = 8, col = "red") + 
    NULL